Publication content

Create datasets

# Load data on individual donations
load("./data/r02.fd.bd.all.rdata") # outputs an object called "output" into the environment
donations <- output

# We only want to look at first donation event values from each donor
donors <- donations %>%
    group_by(donor) %>%
    filter(date == min(date)) %>%
    ungroup()

# Load FinDonor demographic data
load("./data/r02ds.donorData.rdata") # outputs an object called "output" into the environment
findonor <- output

# Combine the FinDonor datasets
FinDonor <- left_join(donors, findonor, by = "donor")

# Load THL data
# Sofie: thldalta.rdata contains all five THL cohorts, extract FINRISK97 and Health2000 from the others
load("./data/thldata.rdata")
FinRisk97 <- thldata$fr1997
Health2k <- thldata$h2000

# Remove leftovers
rm(output)
rm(thldata)

## Rename useful stuff
# Ferritin, Self-Reported Health
FinDonor <- rename(FinDonor, SRH = QR17, Menstruation = QR79, Age_float = Age, Age = age)
FinRisk97 <- rename(FinRisk97, Ferritin = FERRITIN, SRH = Q40, Gender = SUKUP, Menstruation = K129, Age = IKA)
Health2k <- rename(Health2k, Ferritin = FERRITIINI, SRH = BA01, Gender = SP2, Menstruation = BD03, Age = IKA2, Menopause = MENOP, APOB = NMR_APOB, APOA1 = NMR_APOA1)

# Make "useful stuff" conform with each other
FinDonor <- FinDonor %>% 
    mutate(SRH = case_when(SRH == "Excellent" ~ 1,
                           SRH == "Very_good" ~ 2,
                           SRH == "Good" ~ 3,
                           SRH == "Satisfactory" ~ 4,
                           SRH == "Poor" ~ 5),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == "regular_period" | Menstruation == "irregular_period") ~ "Women|Pre",
                             Gender == "Women" & Menstruation == "no_period" ~ "Women|Post",
                             TRUE ~ "NA")) # Equates to "else"

FinRisk97 <- FinRisk97 %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Pre",
                             Gender == "Women" & Menstruation == 3 ~ "Women|Post",
                             TRUE ~ "NA"))

Health2k <- Health2k %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Pre",
                             Gender == "Women" & (Menstruation == 3 | Age >= 55) ~ "Women|Post",
                             TRUE ~ "NA"))

# Donation eligibility
# These are both "approximates" in a sense, we don't have all the necessary variables to
# filter thoroughly, and we'll be able to do more filtering on Health2000 than FinRisk97
donor_eligible_h2k <- Health2k %>%
    filter(BMII_PAINO.x >= 50 | BMII_PAINO.x <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 | Age <= 66) %>% # Filter away too young and too old
    filter((B_Hb >= 125 & Gender == "Women") | (B_Hb >= 135 & Gender == "Men")) %>% # Filter by hemoglobin
    filter(BA08 == 0) %>% # filter out people with heart attacks
    filter(BA09 == 0) %>% # filter out people with angina
    filter(BA10 == 0) %>% # cardiac insufficiency / heart failure
    filter(!(BA26 == 1 & ATC_A10A == 1)) %>% # filter out people who are diabetic AND use insulin
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) %>% # rename for ease of use
    mutate(HbA1C = B_GHb_A1C * 10.93 - 23.50)

donor_eligible_fr <- FinRisk97 %>%
    filter(PAINO >= 50 | PAINO <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 | Age <= 66) %>% # Filter away too young and too old
    #filter((HGB >= 125 & Gender == 2) | (HGB >= 135 & Gender == 1)) %>% # DON'T filter by hemoglobin, < 500 values in data
    filter(Q15A != 2) %>% # STEMI, NSTEMI
    filter(Q16A != 2) %>% # Stroke
    # filter(Q38 != 2 & Q38 != 4) %>%  # Insulin treatment (2: just insulin, 4: insulin and a tablet)
    filter(Q17B != 2) %>% # cardiac insufficiency
    filter(Q17C != 2) %>% # angina pectoris
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) # rename for ease of use

# Create useful mastersets
# fer_srh <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "SRH", "Group")], 
#                      Health2k = donor_eligible_h2k[, c("Ferritin", "SRH", "Group")], .id = "Cohort") %>% 
#     mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
#            Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k")),
#            SRH = ordered(SRH, levels = 1:5)) %>%
#     filter(Group != "NA") %>%
#     drop_na()

fer_crp <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "CRP")], 
                         Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "CRP")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    filter(CRP >= 0.01) %>%
    drop_na()

Mock Cohort|Var Tables

table1 <- as.data.frame(table(fer_crp$Group, fer_crp$Cohort))
table1$CRP <- c(paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1$FER <- c(paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Pre" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Post" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1

##         Var1      Var2 Freq                 CRP                     FER
## 1  Women|Pre FinRisk97 1980 0.76 | (0.38, 1.85)  23.97 | (12.35, 42.51)
## 2 Women|Post FinRisk97 1001  1.26 | (0.61, 2.6)   55.8 | (31.42, 94.74)
## 3        Men FinRisk97 2801 0.91 | (0.47, 1.92) 110.99 | (65.87, 181.3)
## 4  Women|Pre  Health2k  976 0.61 | (0.26, 1.78)   27.9 | (15.17, 48.54)
## 5 Women|Post  Health2k  785  1.02 | (0.4, 2.35)      55.97 | (32, 95.1)
## 6        Men  Health2k 1793 0.79 | (0.35, 1.77)  123.1 | (76.32, 192.7)

Exclusions: Health2000

h2k <- "Health2000"
fr <- "FinRisk1997"
women_m <- "Menstruating women"
women_nm <- "Non-menstruating women"
men <- "Men"
paste0(h2k, ", ", women_m, ": ",  nrow(Health2k %>% filter(Group == "Women|Pre")))

## [1] "Health2000, Menstruating women: 1508"

paste0(h2k, ", ", women_m, ", under 50kg or over 200kg: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(BMII_PAINO.x < 50 | BMII_PAINO.x > 200)))

## [1] "Health2000, Menstruating women, under 50kg or over 200kg: 48"

paste0(h2k, ", ", women_m, ", under 18y or over 66y: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(Age < 18 | Age > 66)))

## [1] "Health2000, Menstruating women, under 18y or over 66y: 0"

paste0(h2k, ", ", women_m, ", Hb under 125: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(B_Hb < 125)))

## [1] "Health2000, Menstruating women, Hb under 125: 284"

paste0(h2k, ", ", women_m, ", heart attack: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(BA08 == 1)))

## [1] "Health2000, Menstruating women, heart attack: 1"

paste0(h2k, ", ", women_m, ", angina: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(BA09 == 1)))

## [1] "Health2000, Menstruating women, angina: 2"

paste0(h2k, ", ", women_m, ", heart failure/insuff: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(BA10 == 1)))

## [1] "Health2000, Menstruating women, heart failure/insuff: 3"

paste0(h2k, ", ", women_m, ", diabetics on insulin: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter((BA26 == 1 & ATC_A10A == 1))))

## [1] "Health2000, Menstruating women, diabetics on insulin: 7"

paste0(h2k, ", ", women_m, ", bad or very bad self-reported health: ",  nrow(Health2k %>% filter(Group == "Women|Pre") %>% filter(SRH >= 4)))

## [1] "Health2000, Menstruating women, bad or very bad self-reported health: 51"

paste0(h2k, ", ", women_nm, ": ",  nrow(Health2k %>% filter(Group == "Women|Post")))

## [1] "Health2000, Non-menstruating women: 1505"

paste0(h2k, ", ", women_nm, ", under 50kg or over 200kg: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(BMII_PAINO.x < 50 | BMII_PAINO.x > 200)))

## [1] "Health2000, Non-menstruating women, under 50kg or over 200kg: 27"

paste0(h2k, ", ", women_nm, ", under 18y or over 66y: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(Age < 18 | Age > 66)))

## [1] "Health2000, Non-menstruating women, under 18y or over 66y: 279"

paste0(h2k, ", ", women_nm, ", Hb under 125: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(B_Hb < 125)))

## [1] "Health2000, Non-menstruating women, Hb under 125: 117"

paste0(h2k, ", ", women_nm, ", heart attack: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(BA08 == 1)))

## [1] "Health2000, Non-menstruating women, heart attack: 30"

paste0(h2k, ", ", women_nm, ", angina: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(BA09 == 1)))

## [1] "Health2000, Non-menstruating women, angina: 90"

paste0(h2k, ", ", women_nm, ", heart failure/insuff: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(BA10 == 1)))

## [1] "Health2000, Non-menstruating women, heart failure/insuff: 52"

paste0(h2k, ", ", women_nm, ", diabetics on insulin: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter((BA26 == 1 & ATC_A10A == 1))))

## [1] "Health2000, Non-menstruating women, diabetics on insulin: 19"

paste0(h2k, ", ", women_nm, ", bad or very bad self-reported health: ",  nrow(Health2k %>% filter(Group == "Women|Post") %>% filter(SRH >= 4)))

## [1] "Health2000, Non-menstruating women, bad or very bad self-reported health: 175"

paste0(h2k, ", ", men, ": ",  nrow(Health2k %>% filter(Group == "Men")))

## [1] "Health2000, Men: 2944"

paste0(h2k, ", ", men, ", under 50kg or over 200kg: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(BMII_PAINO.x < 50 | BMII_PAINO.x > 200)))

## [1] "Health2000, Men, under 50kg or over 200kg: 3"

paste0(h2k, ", ", men, ", under 18y or over 66y: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(Age < 18 | Age > 66)))

## [1] "Health2000, Men, under 18y or over 66y: 229"

paste0(h2k, ", ", men, ", Hb under 125: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(B_Hb < 125)))

## [1] "Health2000, Men, Hb under 125: 40"

paste0(h2k, ", ", men, ", heart attack: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(BA08 == 1)))

## [1] "Health2000, Men, heart attack: 116"

paste0(h2k, ", ", men, ", angina: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(BA09 == 1)))

## [1] "Health2000, Men, angina: 144"

paste0(h2k, ", ", men, ", heart failure/insuff: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(BA10 == 1)))

## [1] "Health2000, Men, heart failure/insuff: 67"

paste0(h2k, ", ", men, ", diabetics on insulin: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter((BA26 == 1 & ATC_A10A == 1))))

## [1] "Health2000, Men, diabetics on insulin: 50"

paste0(h2k, ", ", men, ", bad or very bad self-reported health: ",  nrow(Health2k %>% filter(Group == "Men") %>% filter(SRH >= 4)))

## [1] "Health2000, Men, bad or very bad self-reported health: 271"

Exclusions: FinRisk1997

paste0(fr, ", ", women_m, ": ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre")))

## [1] "FinRisk1997, Menstruating women: 2469"

paste0(fr, ", ", women_m, ", under 50kg or over 200kg: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(PAINO < 50 | PAINO > 200)))

## [1] "FinRisk1997, Menstruating women, under 50kg or over 200kg: 81"

paste0(fr, ", ", women_m, ", under 18y or over 66y: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(Age < 18 | Age > 66)))

## [1] "FinRisk1997, Menstruating women, under 18y or over 66y: 8"

paste0(fr, ", ", women_m, ", STEMI, NSTEMI: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(Q15A == 2)))

## [1] "FinRisk1997, Menstruating women, STEMI, NSTEMI: 6"

paste0(fr, ", ", women_m, ", stroke: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(Q16A == 2)))

## [1] "FinRisk1997, Menstruating women, stroke: 11"

paste0(fr, ", ", women_m, ", heart insuff: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(Q17B == 2)))

## [1] "FinRisk1997, Menstruating women, heart insuff: 12"

paste0(fr, ", ", women_m, ", angina: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(Q17C == 2)))

## [1] "FinRisk1997, Menstruating women, angina: 37"

paste0(fr, ", ", women_m, ", bad or very bad self-reported health: ",  nrow(FinRisk97 %>% filter(Group == "Women|Pre") %>% filter(SRH >= 4)))

## [1] "FinRisk1997, Menstruating women, bad or very bad self-reported health: 114"

paste0(fr, ", ", women_nm, ": ",  nrow(FinRisk97 %>% filter(Group == "Women|Post")))

## [1] "FinRisk1997, Non-menstruating women: 1463"

paste0(fr, ", ", women_nm, ", under 50kg or over 200kg: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(PAINO < 50 | PAINO > 200)))

## [1] "FinRisk1997, Non-menstruating women, under 50kg or over 200kg: 31"

paste0(fr, ", ", women_nm, ", under 18y or over 66y: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(Age < 18 | Age > 66)))

## [1] "FinRisk1997, Non-menstruating women, under 18y or over 66y: 156"

paste0(fr, ", ", women_nm, ", STEMI, NSTEMI: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(Q15A == 2)))

## [1] "FinRisk1997, Non-menstruating women, STEMI, NSTEMI: 34"

paste0(fr, ", ", women_nm, ", stroke: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(Q16A == 2)))

## [1] "FinRisk1997, Non-menstruating women, stroke: 43"

paste0(fr, ", ", women_nm, ", heart insuff: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(Q17B == 2)))

## [1] "FinRisk1997, Non-menstruating women, heart insuff: 69"

paste0(fr, ", ", women_nm, ", angina: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(Q17C == 2)))

## [1] "FinRisk1997, Non-menstruating women, angina: 94"

paste0(fr, ", ", women_nm, ", bad or very bad self-reported health: ",  nrow(FinRisk97 %>% filter(Group == "Women|Post") %>% filter(SRH >= 4)))

## [1] "FinRisk1997, Non-menstruating women, bad or very bad self-reported health: 179"

paste0(fr, ", ", men, ": ",  nrow(FinRisk97 %>% filter(Group == "Men")))

## [1] "FinRisk1997, Men: 3943"

paste0(fr, ", ", men, ", under 50kg or over 200kg: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(PAINO < 50 | PAINO > 200)))

## [1] "FinRisk1997, Men, under 50kg or over 200kg: 4"

paste0(fr, ", ", men, ", under 18y or over 66y: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(Age < 18 | Age > 66)))

## [1] "FinRisk1997, Men, under 18y or over 66y: 375"

paste0(fr, ", ", men, ", STEMI, NSTEMI: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(Q15A == 2)))

## [1] "FinRisk1997, Men, STEMI, NSTEMI: 160"

paste0(fr, ", ", men, ", stroke: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(Q16A == 2)))

## [1] "FinRisk1997, Men, stroke: 95"

paste0(fr, ", ", men, ", heart insuff: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(Q17B == 2)))

## [1] "FinRisk1997, Men, heart insuff: 155"

paste0(fr, ", ", men, ", angina: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(Q17C == 2)))

## [1] "FinRisk1997, Men, angina: 224"

paste0(fr, ", ", men, ", bad or very bad self-reported health: ",  nrow(FinRisk97 %>% filter(Group == "Men") %>% filter(SRH >= 4)))

## [1] "FinRisk1997, Men, bad or very bad self-reported health: 382"

Ferritin X CRP

Subgroups: menstruating women, non-menstruating women, men

options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    theme_minimal() + 
    geom_smooth(method = "lm", formula = y ~ x, color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort))

Ferritin X CRP | COLORED

options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

Ferritin X CRP | GROUPS ONLY BY COLOR

options(scipen = 10000)
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    geom_point(aes(color = Group), alpha = 0.2) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(cols = vars(Cohort)) +
    theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(alpha = 1)))

Proportion analysis code

The estimate: How many more people are over 3 mg/l in CRP when we filter the populations by different ferritin thresholds? Expressed in percentage point difference. The intuition: We currently don’t do ferritin filtering. The population from which we draw donors has a certain proportion of people in risk of cardiac/chronic inflammation, often indicated by elevated hs-CRP (>3 mg/l). Does this proportion increase in a statistically significant manner, if we start requiring higher ferritin levels from donors?

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
CRP_trld <- 3

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_crp %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_crp %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i],var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_CRP", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))) { # run bootstrap only if needed
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_crp %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_crp %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_CRP", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))

Proportion analysis plot

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + guides(linetype = "none")

Checking statistical significance at points of interest

IMPORTANT: We are checking for non-overlaps of 95% confidence intervals. These intervals have been derived from 100 bootstrap samples using normal approximation. Normal approximation relies on the assumption, that our bootstrapped estimates are normally distributed, which we have confirmed separately using both visual assessment (a histogram) and a standard Shapiro-Wilk test.

# FINRISK
# Menstruating Women
frwomenpre5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 5)
frwomenpre15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 15)
frwomenpre30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 30)
frwomenpre50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Pre" & Ferritin == 50)

# Non-menstruating Women
frwomenpost5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 5)
frwomenpost15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 15)
frwomenpost30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 30)
frwomenpost50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Post" & Ferritin == 50)

# Men
frmen5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 5)
frmen15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 15)
frmen30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 30)
frmen50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 50)

# H2K
# Menstruating Women
h2kwomenpre5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 5)
h2kwomenpre15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 15)
h2kwomenpre30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 30)
h2kwomenpre50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Pre" & Ferritin == 50)

# Non-menstruating Women
h2kwomenpost5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 5)
h2kwomenpost15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 15)
h2kwomenpost30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 30)
h2kwomenpost50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Post" & Ferritin == 50)

# Men
h2kmen5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 5)
h2kmen15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 15)
h2kmen30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 30)
h2kmen50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 50)

Using these objects, we can check the significance of between the points of interest. For example “h2kwomenpre5\(upper < h2kwomenpre15\)lower” would evaluate to TRUE, if the proportion of menstruating women over 3 mg/l CRP is significantly higher in the population filtered by ferritin 15 ug/l than 5 ug/l. Indeed, for menstruating women, the proportion is significantly higher at filter levels of 15, 30, and 50 ug/l when compared with 5 ug/l. This holds for non-menstruating women also, except for the FinRisk 1997 cohort, where the difference between levels 5 and 15 was not significant. The differences are all significant also in men, but the respective increases in proportions are much smaller.

Proportion analysis plot | COLORED

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") +
    theme(legend.position = "none")

Proportion analysis plot | COLORED & ALL-IN-1

ggplot(data = means_all, aes(x = Ferritin, y = means, group = Group)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group, shape = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(cols = vars(Cohort)) +
    labs(y = "%p") + 
    guides(fill = "none") +
    theme(legend.position = "bottom")

Additional markers (for the supplement)

We’ll also take a look at acetylated glycoprotein measurements, glucated haemoglobin, and apolipoproteins A1 and B.

# mastersets for the Supplement
# GlycA
fer_glyca <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "GlycA")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "GlycA")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# HbA1C
fer_hba1c <- bind_rows(Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "HbA1C")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()

# APOB
fer_apob <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOB")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOB")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# APOA1
fer_apoa1 <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOA1")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOA1")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Pre", "Women|Post", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()

Mock Cohort|Var Tables (Supplement)

suptable1 <- as.data.frame(table(fer_glyca$Group, fer_glyca$Cohort))

suptable1$GlycA <- c(paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable1$FER <- c(paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Pre" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Post" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable2 <- as.data.frame(table(fer_hba1c$Group))
suptable2$HbA1C <- c(paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable2$FER <-  c(paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Pre" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Post" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable1

##         Var1      Var2 Freq               GlycA                     FER
## 1  Women|Pre FinRisk97 2045 1.28 | (1.17, 1.42)     23.9 | (12.3, 42.5)
## 2 Women|Post FinRisk97 1002 1.38 | (1.26, 1.53)  55.79 | (31.42, 94.53)
## 3        Men FinRisk97 2821 1.39 | (1.26, 1.56) 111.37 | (66.46, 181.9)
## 4  Women|Pre  Health2k 1108 1.09 | (0.95, 1.23)    27.6 | (14.8, 48.23)
## 5 Women|Post  Health2k  834     1.15 | (1, 1.3)   55.8 | (31.42, 94.36)
## 6        Men  Health2k 1990 1.19 | (1.05, 1.35) 120.8 | (75.15, 189.15)

suptable2

##         Var1 Freq                  HbA1C                   FER
## 1  Women|Pre 1111 31.15 | (28.96, 33.34)  27.6 | (14.8, 48.17)
## 2 Women|Post  836 33.34 | (31.15, 36.61) 55.8 | (31.48, 94.51)
## 3        Men 1997 34.43 | (32.24, 36.61) 120.8 | (75.5, 189.2)

# TODO: suptable3
# TODO: suptable4

Ferritin X GlycA | COLORED

options(scipen = 10000)
ggplot(data = fer_glyca, aes(x = Ferritin, y = GlycA)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

Proportion analysis code (GlycA)

For the acetylated glycoprotein we don’t have any well established thresholds for “healthy” and “unhealthy”. We’ll use the population median here.

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "FinRisk97", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_glyca %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_glyca %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_glyca", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "Health2k", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_glyca %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_glyca %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_glyca", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))

Proportion analysis plot (GlycA)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")

Ferritin X HbA1C | COLORED

options(scipen = 10000)
ggplot(data = fer_hba1c, aes(x = Ferritin, y = HbA1C)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group)) +
    theme(legend.position = "none") +
    labs(y = expression(HbA[1*C]))

Proportion analysis code (HbA1C)

The reference values for healthy people are between 20 - 42 mmol/mol (https://www.terveyskirjasto.fi/snk03092). Because the relationship with ferritin appears to be positive, we’ll use the upper bound.

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
HbA1C_trld <- 42

if (!file.exists(paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_hba1c %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_hba1c %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_hba1c %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_hba1c", boot_n, ".rds"))}

means_all <- means_health2k
means_all$Cohort <- c(rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))

Proportion analysis plot (HbA1C)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group)) +
    labs(y = "%p") + theme(legend.position = "none")

Ferritin X APOB | COLORED

options(scipen = 10000)
ggplot(data = fer_apob, aes(x = Ferritin, y = APOB)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

Ferritin X APOA1 | COLORED

options(scipen = 10000)
ggplot(data = fer_apoa1, aes(x = Ferritin, y = APOA1)) + 
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none")

Proportion analysis code (APOB)

Apolipoprotein B elevates with the risk of cardiovascular diseases and it correlates positively with ferritin. The reference values for healthy men are 0.6 - 1.5 (of which we’ll use 1.5) and for healthy women they are 0.6 - 1.3 (we’ll use 1.3). (https://huslab.fi/ohjekirja/20705.html)

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
APOB_trld_men <- 1.5 # for males [0.6, 1.5]
APOB_trld_women <- 1.3 # for females [0.6, 1.3]

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_apob %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_apob %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_APOB", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_apob %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_apob %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_APOB", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))

Proportion analysis plot (APOB)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")

Proportion analysis code (APOA1)

Apolipoprotein A1 decreases with the risk of atherosclerosis (peripheral, coronary, brain) and appears to correlate negatively with ferritin. Thus, we will use the lower bounds of the reference intervals for healthy individuals here. These are 1.1 - 2.0 for men and 1.2 - 2.3 for women. (https://huslab.fi/ohjekirja/20705.html)

ferritin_values <- seq(5, 50, 1)
iterations <- length(ferritin_values)
APOA1_trld_men <- 1.1 # for males [1.1, 2.0]
APOA1_trld_women <- 1.2 # for females [1.2, 2.3]

if (!file.exists(paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        boot_obj_men <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
        ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
        # Women|Pre
        boot_obj_women_pre <- boot(fer_apoa1 %>% filter(Group == "Women|Pre" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
        # Women|Post
        boot_obj_women_post <- boot(fer_apoa1 %>% filter(Group == "Women|Post" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
        
        ## Store
        # Men
        means_men[i] <- boot_obj_men$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Pre
        means_women_pre[i] <- boot_obj_women_pre$t0
        upper_women_pre[i] <- ci_obj_women_pre$normal[3]
        lower_women_pre[i] <- ci_obj_women_pre$normal[2]
        # Women|Post
        means_women_post[i] <- boot_obj_women_post$t0
        upper_women_post[i] <- ci_obj_women_post$normal[3]
        lower_women_post[i] <- ci_obj_women_post$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_pre, means_women_post),
                                     upper = c(upper_men, upper_women_pre, upper_women_post),
                                     lower = c(lower_men, lower_women_pre, lower_women_post),
                                     Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))
    
    }
    
    # Save
    saveRDS(means_finrisk, paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))
} else {means_finrisk <- readRDS(paste0("./data/PUBL_finrisk_ratio_APOA1", boot_n, ".rds"))}

if (!file.exists(paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Pre
    means_women_pre <- 1:iterations
    upper_women_pre <- 1:iterations
    lower_women_pre <- 1:iterations
    # Women|Post
    means_women_post <- 1:iterations
    upper_women_post <- 1:iterations
    lower_women_post <- 1:iterations

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    boot_obj_men <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
    ci_obj_men <- boot.ci(boot_obj_men, type = "norm")
    # Women|Pre
    boot_obj_women_pre <- boot(fer_apoa1 %>% filter(Group == "Women|Pre" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_pre <- boot.ci(boot_obj_women_pre, type = "norm")
    # Women|Post
    boot_obj_women_post <- boot(fer_apoa1 %>% filter(Group == "Women|Post" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_post <- boot.ci(boot_obj_women_post, type = "norm")
    
    ## Store
    # Men
    means_men[i] <- boot_obj_men$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Pre
    means_women_pre[i] <- boot_obj_women_pre$t0
    upper_women_pre[i] <- ci_obj_women_pre$normal[3]
    lower_women_pre[i] <- ci_obj_women_pre$normal[2]
    # Women|Post
    means_women_post[i] <- boot_obj_women_post$t0
    upper_women_post[i] <- ci_obj_women_post$normal[3]
    lower_women_post[i] <- ci_obj_women_post$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_pre, means_women_post),
                                 upper = c(upper_men, upper_women_pre, upper_women_post),
                                 lower = c(lower_men, lower_women_pre, lower_women_post),
                                 Gender = c(rep("Men", iterations), rep("Women|Pre", iterations), rep("Women|Post", iterations)))    
    
    
    }
    
    # Save
    saveRDS(means_health2k, paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))
} else {means_health2k <- readRDS(paste0("./data/PUBL_health2k_ratio_APOA1", boot_n, ".rds"))}

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 138), rep("Health2k", 138))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Pre", "Women|Post", "Men"))

Proportion analysis plot (APOA1)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Post", "Women|Pre" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(y = "%p") + theme(legend.position = "none")